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1.0 INTRODUCTION 


The User's Manual is meant to serve as a guide in helping the user 
make successful runs with the MINT computer program. The degree of success 
obtained by the user will depend on the skill of the user and his ability 
to correctly apply the code to his particular problem. The code will solve 
the governing equations, subject to the boundary conditions. However, 
meaningful results will only be obtained if the boundary conditions are 
appropriate to the problem. In addition, the user must specify viscosity 
models, initial conditions, a coordinate system and the location of the 
grid pointe to adequately resolve the flow. The user with a good knowledge 
of the physics and a moderate amount of experience should be able to 
successfully apply the code to a wide variety of problems. 

The User's Manual is divided into eleven parts consieting of: 
(1) a flow diagram of che MINT compiter code, (2) a brief description of 
each subroutine and fits use, (3) a Liet of the major PORTRAN vartatles and 
a description of their meaning, (4) a description of the logical file units 
utilized by the MINT computer code, (5) a brief discourse on the general 
strategy for the running of the MINT computer code, (6) a description of 
the method used to distribute grid points on the boundaries, (7) the 
protocol for grid point and boundary identification, (8) a detatled 
description of the fnput required to run labyrinth seal cases, (9) sample 
{uput for a labyrinth seal case, (10) sample output for the corresponding 
case and (11) a deseri{ption and condensed aser's manual for a plot prograu 
that can be used with the MINT computer code. For detaile of the equations 


solved, boundary conditions utilized and the numer.cal procedure, the user 


{s referred to Volume tL of this cepord, Ref. 1. 


2.0 PLOW DIAGRAM 


The purpose of the flow diagram is to help the user understand the 
basic flow of information within the MINT computer code. Because of the 
size of the code (approximately 20,000 cards), a detailed flow diagram 
would be prohibitively large and probably of little value to the ':3er. 
Therefore, the flow diagram is intended te give the user a genersl overview 
of the structure of the code. The interested user is urged to consult the 


program listing for details. The subroutines contain comment cards to aid 


the user in the interpretation of the coding. 


START 


Read Input and 
Initialize Flow Field 


Call EXEC 


March Calculation 
; NT Time Steps 


Call RESTRT 


Write Restart Files 


Call PRNTECPRNTA) 


Write Final Flow Field 
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Se eS a a aS foe Sane! 


Called from DAL 


Set Default Input 


Call RDIND1(READ1) 


= <a> 


Yes 


Read NAMELIST 
Input from 
File INTAP1 


Read NAMELIST 
Input from 
File 5 


Set Integer 
Variables and Flags 


Call BVIV 


Set Grid and 
Initialize Flow Field 


Yes Set Grid Parametecs 


No . Call COORD 


ee ee 


Call RDIND2(READB) Set Geometr 
READ Flow Field from 
Restart File INTAP 


Return to DAL 


No 


Yes 


Initialize Flow Field 
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Called from DAL 


Call EXTBV 


ee ee 


Calculate 
DtD V-¥, Viscosity 
Zero Analysis 
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Call EXTBV 
Set Boundary Valu2s 
Update Pressure 

Temperature 
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Call PRNTS 
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IDUMP1=1 


Intermediate Output 


Adjust Time Step 


-—=— es =— me om om 


Flow Field 


Do 
1000"I1 DT 


Return to DAL 


Call TIMGEO 


Set Geometry 


Call ADIC 
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Solve Equations 
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Called from 
ADIC 
Gall ADICX 


at 


Call ADICP1 
Generate Equations 
Obtain Solution 
First ADI Sweep 


Call ADICZI 


Call ADICP2 


~~ me meme 


Generate Equations 
Obtain Solution 
Second ADI Sweep 


Return to ADIC 
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ia 3.0 LIST OF MAJOR FORTRAN VARIABLES 
18 
7 oo FORTRAN COMMON 
. SYMBOL BLOCK DESCRIPTION 
ee AC(I,J) BLK] DEPENDENT VARIABLE ARRAY 
. _ Acs(T,2) BLKi GEOMETRY DATA ARRAY 
i ae ANCL, J) BLKM ARRAY STORING TIME TERM LINEARIZED COEFFICIENTS 
Ny APR(I, J) PRNT PRINT OUTPUT ARRAY 
C - AVANDR TURB DAMPING CONSTANT 
: “ 5 AVISC(1,J) . MISC? ARTIFICIAL DISSIPATION PARAMETER 
: C(1,J,K) BLEM COUPLED MATRIX ARRAY STORAGE 
_ CLENG CREF REFERENCE LENGTH 
CMACH MISC2 REFERENCE MACH NUMBER 
D VARNO INDEX FOR DIVERGENCE 
yO DI(I,J,K) BLE RRA. STORING PIRST SWEEP LINEARIZED COEFFICIENTS 
"a 02(1,J,K) - BLKK ARRAY STORING SECOND SWEEP LINEARLZED COEFFICIENTS 
D3(1,3,K) BLE ..RRAY STORING TsIRD SWEEP LINRARIZED COEFFICIENTS 
. DENSR CREF REF SRENCE DENSITY 
; DFW(I,J,K) ADI? ONFPERENCE WEIGHT ARR/Y 
‘ oy DIMI NOND INVERSE REYNOLDS NUMBFR 
a DIM2 NOND REFERENCE PRESSURL/REPERENCE GYNAMIC HEAD 
4 DIN3 NOND REFERENCE PRESSURE/(REPERENCE DENSITY * 
: REFGRENCE ENTHALPY) 
x DIMA NOND 1.0/(REY * PRNCL) 
. DIM12 NOND 2.0 * DIMI 
" : DS VARNO so INDEX FOR DISSIPATION 
‘ ot MISC2 «TIME STEP 
7 
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FORTRAN 
SYMBOL 


DTCON 
DTMAX 
DIMIN 
E(I,J,K) 
GAMMA 

H 

Tl 

ITADI 


IBC 


IBOUND(I, JDIR) 


ICOORD 
IDT 
IDTADJ 


IDUMP1 


1EQ 
LEQBC(1,J,K) 
1FSB 
IGPRT(I) 
IL 

LNOUT 
IPRINT 
IREST 
IUNITS 
IVARPR( 1) 
JADI 


JG 


COMMON 
BLOCK 


MISC2 


MISC2 


MISC2 


BLEM 


CREF 


VARNO 


MGAUS 


ADI1 


ADI1 


BCTYP 


GEO] 


MISC2 


MISC2 


OUTA 


ADI1 


ADI} 


PRNT 


GEO! 


NGAUS 


CMDISK 


MISC2 


MISC2 


CNEF 


NISC2 


ADI1 


None 


DESCRIPTION 

INVERSE STEP 

MAXIMUM ALLOWABLE TIME STEP (SEE $READY) 
MINIMUM ALLOWABLE TIME STEP (SEE $READY) 
COUPLED MATRIX ARRAY STORAGE 

RATIO OF SPECIFIC HEATS 

INDEX FOR ENTHALPY 

LOWER LIMIT FOR MATRIX INVERSION 

ADI SWEEP NUMBER 

BOUNDARY CONDITION BOUNDARY PARAMETER 
BOUNDARY TYPE 

COORDINATE SYSTEM JNDEX (SEE $READ2) 
TIME STEP INDEX 

TIME STEP CONTROL PARAMETER (SEE $READ4) 


PARAMETER CONTROLLING INITIAL STATION 
PRINT 


EQUATION NUMBER 

BOUNDARY CONDITION INDICATION (SBE $READ4) 
GRID INDICATOR 

GEOMETRY PRINT CONTROL (SEE $READ4) 

UPPER LIMIT FOR MATRIX INVERSION 

IN CORE - OUT OF CORE OPTION 

PRINT INTERVAL PARAMETER (SEE $READ4) 
RESTART READ CONTROL PARAMETER (SEE SREAD1) 
DIMENSIONAL UNITS SENTINEL (SEE SREAD3) 
PRINT PARAMETER (SPE S$READ4) 

ADI SWEEP PARAMETER 


IMPLICIT GRID POINT NUMBER 


FORTRAN 
SYMBOL 


JX 
KZ 
LCORNR 
LCORNX (1) 
LCORNZ (1) 


LX 


REY 
SAVE 


COMMON 
BLOCK 


ADI2 


ADI2 


DBLVAL 


DBLVAL 


DBLVAL 


ADI2 


ADI} 


- MISC2 


REDUMP 
MISC2 
MISC2 
MISC2 
ADI4 
ADIG 
ADIG 
ADI4 
VARKO 
MISC2 
MISC2 
CREF 
CREF 


CREF 


DESCRIPTION 

DIRECTION-1 GRID POINT INDEX 

DIRECTION-3 GRID POINT INDEX 

NUMBER OF RE-ENTRANT CORNERS 

X-LOCATION OF RE-ENTRANT CORNER 

Z-LOCATION OF RE-ENTRANT CORNER 

DIRECTION-! GRID POINT INDEX ~ 

NUMBER OF EQUATIONS TO BE SOLVED 

NUMBER OF TIME STEPS TO BE RUN 

NUMBER OF TIME STEPS BETWEEN RESTARTS (SEE $READ4) 
NUMBER OF TIME STEPS IN A CYCLE (SEE $READ4) 
NUMBER OF INTERIOR DIRECTION-1 POINTS (SEE $READ2) 
NUMBER OF INTERIOR DIRECTION-3 POINTS (SEE $READ2) 
FIRST GRID POINT - DIRECTION 1 

LAST GRID POINT - DIRECTION 1 

FIRST GRID - DIRECTION 3 

LAST GRID POINT ~ DIRECTION 3 

INDEX FOR PRESSURE 

TIME STEP CONTROL PARAMETER 

TIME STEP CONTROL PARAMETER 

DIMENSIONAL FREE STREAM PRESSURE (SEE SREAD3) 
REFERENCE PRESSURE 

PRANDIL NUMBER 

TOTAL PRESSURE 

DIMENSIONAL STAGNATION PRESSURE (SEE $READ3) 

INDEX FOR DENSITY 

REYNOLDS NUMBER PER UNIT LENGTH (SEE SREAD3) 


REYNOLDS NUMBER 
STORAGE FOR * LEVEL DELTA VALUES 
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FOTRAN COMMON 
SYMBOL BLOCK DESCRIPTION 
SN(I) BLKM ARRAY STORING SOURCE TERM LINEARIZED COEFFICIENT 
SSTEST MISC2 MAXIMUM CHANGE IN VARIABLE ACROSS TIME STEP 
T VARNO —_—sCINDEX FOR TEMPERATURE 
TAUW TURB WALL SHEAR 
TED VARNO —sINDEX FOR TURBULENCE KINECTIC ENERGY 
TKE VARNO —sCINDEX FOR DISSIPATION OF TURBULENCE KINECTIC ENERGY 
TREF - CREF REFERENCE TEMPERATURE 

, TTIME MISC2 CUMULATIVE TIME 

x TTOT BCCON TOTAL TEMPERATURE 

‘ U VARNO —sCINDEX FOR DIRECTION~1 VELOCITY 

; USTAR TURB DIMENSIONLESS VELOCITY 

v VARNO —sINDEX FOR DIRECTION~2 VELOCITY 

: 4 VISCL TURB LAMINAR REFERENCE VISCOSITY 

a : VISCR CREF REFERENCE VISCOSITY 

; vS1 VARNO —s INDEX FOR LAMINAR VISCOSITY 

f . vs2 VARNO —s INDEX FOR EFFECTIVE VISCOSITY 

\ : w VARNO  —sINDEX FOR DIRECTION-3 VELOCITY 

: , WREF CREF REFERENCE VELOCITY 

5 § XGMAX (1) GRIDI MAXIMUM COORDINATE VALUE (SEE $READ2) 

; : XCHIN(T) GRID MINIMUM COORDINATE VALUE (SEE $READ2) 

td 4 YPLUS TURB DIMENSIONLESS DISTANCE FROM SURFACE 

mn 

. 

t ‘ 

X 
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Subroutine 


ADIC 


ADICPI 


ADICP2 


ADICX 
ADICZ1 
ADICZO 


ADIUN1 


ADIUN2 


AMATRX 
ARTVIS1 
ATIME 


BC 


BCPN 


WIV 


CONVCT 1 


CONVCT2 


COORD 
CORNER 
CURVT 
DAL 


DELU 


4.0 MINT SUBROUTINES 
Purpose 
Primary control subroutine for ADI procedure. 


Control subroutine for first ADI sweep solution - coupled 
equations. 


Control subroutine for last ADI sweep solution - coupled 
equations. 


Control subroutine for first ADI sweep. 
Control subroutine for last ADI sweep — in core option. 
Control subroutine for last ADI sweep — out of core option. 


Control subroutine for first ADI sweep solution ~ uncoupled 
equations. 


Control subroutine for last ADI sweep solution — uncoupled 
equations. 


Control eubroutine for time and source terms. 
Calculates artificial dissipation terns. 
Linearizes time derivatives. 


Boundary condition subroutine for general boundary 
conditions. 


Spectalized boundary conditions cubroutine for pressure and 
soacntua. 


Initiates finite difference grid and geometric parameters. 


Lineariees convection teras for cartesian velocity component 
formulation. 


Linearizes convection teras for cylindrical polar velocity 
component formulation. 


Initializes geometry. 

Sets up logic for rce~entrant corners. 
Linearlzes cutvature terms. 

Main prograa. 


Calculates dilatation tetas. 
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Subrout ine 


DIFF 1 


DIFF2 


DISFCN 


EOSDP 


EOSUPI 


EXEC 


EXTBY 


FLWFLD 
GASP 
GAUSS 


GENBC 


GENEQI 


GENEQ?2 


GENUBC 


GEOQN 
GSUPS 


INTGEO 


Purpose 


Linearizes diffusion terms for cartesian velocity component 
formulation. 


Linearizee diffusion terms for cylindrical polar velocity 
component formulation. 


Calculates dissipation function. 
Calculates divergence of the velocity. 
Control subroutine for partial differential equation writing. 
Determines logic for bookkeeping of re~entrant corner logic. 
General equation of state subroutine. 

Linearizes pressure gradient terms. 


Updates pressure and temperature; calculates density and 
enthalpy given pressure and temperature. 


Master control subroutine for time marching of solution. 
Subroutine which updates solution variables from the nt® co 
the n¢18° time level and calculates corresponding values 
for pressure, Cemperature, viscosity, etc. 

Subroutine to generate initial guess. 

Subroutine which calculates properties of a gas. 


Solves the tridiagonal matrix equation Ax@B for x. 


Control subroutine for boundary conditions - coupled 
equations. 


Partial differential equation writing subroutine for 
cartesian velocity component formulation. 


Partial differential equation writing subroutine for 
cylindrical polar velocity component formulation. 


Control subroutine for boundary conditions ~ uncoupled 
equations. 


Geometric initialization subroutine. 
Calculates nonorthogonal traneformation information, 


Geometcic initialization subroutine. 
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he Subrout ine Purpose 
: LOADAC Out of core control subroutine. 
z 
: MGAUSF Block tridiagonal matrix solver. 
% 
: MINMAX Calculates minimum and maximum extents of a domain. 
iF ; . : 
4 MINV Subroutine which solves general matrix equation Ax=B for x. 
R MIXLEN Calculates mixing length. 
hb 
ie NORMD Normalization subroutine. 
; OUGR ID Distributes grid points on a line by the use of error 
. * functions - method of Oh. 
la PRGEO Controls printout of geometric variables. 
i PRINT} Prints out fluid dynamic variables. 
is PRNTA Master control subroutine for printout of fluid dynamic 
te ; variables. 
b ’ 
3 PROB Calculates probability distribuiton function. 
fe 
rm QUICK Quick matrix solver. 
nd 
i RDLIST Reade and writes NAMELIST input data. 
* 
Ma READA Master control eubroutine reads in input and initializes 
flow. 
x READB Reads in data and positions unite for restart files. 
K 
y RESTRT Controle reading and writing of restart files. 
RESULT Control subroutine for printer plots. 
ROTATE Shifts iaformation for out of core option. 
a SETBV Control subroutine for performing endcap boundary conditions. 
oa SETDKE Calculates initial condition for k-e turbulence model. 
ee 
S sorcil Calculates for Reynolds number correction terus-cartesian 
a foraul at ior. 
& 
he ; ‘ 
ee SORCH2 Calculates for Reynolds number correction terms - cylindrical 
Eq polar formulation. 
‘ SOURCE General subroutine which calculates source terms. 
Ls ; 
iN SSTST Calculates maximum charges during a Cime step, 
i 
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Subroutine 


STREAM 
TANHYP 
TEMPN 


TIMGEO 


TIMLFT 


TPLOT 


TURBP 


UPDRUM 


VISCOS 


VKPOL 


WALLFN 


WHERE 


WRDISK 


WRPLOT 


WRSLAB 


Purpose 


Calculates stream function. 
Hyperbolic tangent distribution of grid points along a line. 
Control subroutine to update temperature and pressure. 


Primary subroutine which calculates geometric transformation 
information. 


Subroutine which calculates CPU time remaining. 

Calculates information needed for printer time plots. 
Calculates turbulent profile by method of Maise and McDonald. 
Disk control subroutine. 

Calculates laminar and turbulent viscosity. 

Von Karman - Polhausen velocity profile generator. 

Calculates velocity profile information for rough walls. 
Subroutine which calculates location of nearest walls. 


General purpose diak file read, write and eee 
subroutine. 


Subroutine which writes plot file. 


General putpose disk manipulation subroutine. 
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5.0 LOGICAL FILE UNITS UTILIZED BY THE MINT COMPUTER CODE 


The MINT computer code utilizes up to ten (10) logical file units 
during the execution of a run. In many cases not all ten (10) units are 
used, and hence in these cases it {s not necessary to define all ten (10) 
units. All references to logical file units in the MINT computer code are 
accomplished through the use of a FORTRAN name rather than through a 
specific unit number. Thus if the user desires to change a logical file 
unit number, this can be accomplished easily. A list of the logical file 
units utilized by the MINT computer code, their FORTRAN names, the default 
value unit number and a brief description of the use of the unit is 


presented below. 


FORTRAN Name Default Unit Number Description 

KTAPE 1 Plot file. 

MINP 5 Input file for NAMELIST Input. 
MPRT 6 Printed output unit. 

IOTAPS 10 Output unit nusber for 


dependent variable array. 


INTAPE 10 Input unit number for 
dependent variable array. 


NUNITI ul Scratch file for out’ >f 
core option. 


NUNIT2 12 Seratch file for out of 
core option. 


IOTAP) 20 Output unit number for 
NAMELIST reetart data. 


INTAP! 20 Input unit number for 
NAMELIST restart data. 


NTAG 21 Input unit for geometric 
data - X,Z pairs vritter 
in records of length NDUMX+42. 
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6-0 STRATEGY FOR RUNNING THE MINT COMPUTER CODE 


The following discussion is meant to give the user a basic 
understanding of the strategy for successfully running the MINT computer 
code. It is suggested that the user read and understand this section 
before attempting to prepare input for the MINT code. Basically the 
strategy for running the MINT code can be divided into seven areas: 

(1) Specification of the system of unite to be employed, and reference flow 
conditions, (2) generation of the coordinate system in which the 
calculation is to be performed, (3) specification of the boundary 
conditions to be used, (4) generation of the initial conditions 

(5) selection of the time step to insure rapid convergence of the solution 
(6) solution of the governing equations with a k-€ turbulence model and 

(7) solution of the governing equations with rotation. 

Preparation of the input for the MINT computer code first requires 
the choice of a system of units to be used for the nondimensionalization of 
the governing equations (sll variabies are nondimensionalized with respect 
to the free stream conditions). Two options are available in the MINT 
code, viz., the MKS eystem and the English system corresponding to values 
of the FORTRAN variable IUNITS« 1 and 2 respectively. Table 1 shows the 
various reference variables, and the unita in both the MKS and English 
syetems. Reference condition input for the MINT code can be divided into 
two categories which are referred to av: (1) flight conditions and (2) wind 
tunnel conditions and tefer to the types of Flow condition inforaation that 
would noreaally be available to the user. Adequate discusston of the two 
types of reference conditions requires the use of several simple 
relationships which tle together the various flow conditions. They are the 
perfect gas law 


Py = pPeRT (1) 
the definition of Mach auaber 
Uo 
rey (2) 
- JyRT, 


the isentropic relationships along a streamline for a perfect gas with a 


constant specific heat, Cyypo » 
(3) 
yr 


Pp 
gan. Jo! 2 


(4) 
TO ay eget gs 
Te 252 


The definition of Reynolds number per unit length 


Re. Peace 


pa Di en 5 
t., BeolTeo) 


Sutherlands viscosity law, which asaumes that the viscosity of a 


thermally perfect gas 1s a function of temperature (Te) alone, 


eye 
Hea (Teo) = ~ ” 
G+ T, 


the specific heat relationship for a perfect gaa 
Y 
fe. foe (7) 
end the definition of enthalpy for a conatant epectfic heat perfect gas, 


hy, = op Te (8) 

When flight conditions are given, typically a reference length, Le, 
free stream Mach nuaber, Me, static tesperatute, Te, and static 
ptessutre, Pe, ate specified. Given these variables the free e* ream or 
reference velocity, Us, can be calculated from Bq. (2), while Eq. (1) can 
be used to calculate the free etreas density, Oe. Eq. (6) yields the 
free stream viscosity, Us, and Eq. (5) can be used to calculate the 
Reynolds number per unit liscgth. qe. (3) end (4) yleld the free etreas 
etagnation pressure and tomperature while Eqs. (7) and (8) yteld the free 


strean specific heat and enthalpy. 
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When wind tunnel conditions are given, the normal input information 
usually consists of the refrence length, oo, Mach vumber, Mo, the 
Reynolds number per unit length, Re/Lo, and ef.uer the static pressure, 
Po, or the stagnation pressure, Py. By combining Eqs. (1), (2), (5) 
and (6) a quadratic equation for Tw cea »e obtained in terms of known 


variables, i.e., 


1. Bta/BPt deyeqi2 (9) 


ie 2c, 
where 
Rh 
sive R M Po . R MP) 
~ Somme ry Looe 
Re Re (1+ 22!y2)" y-1 (10) 
lw i 2% 


Eqs. (1), (2), (4), (6), (7) and (8) can then be used to calculate free 
stream density, velocity, stagnation temperature, viscosity, specific heat 
and enthalpy respectively. 

Other categories of input are also possible and any of the free 
stream or reference conditious can be calculated by the appropriate 
combination of Eqs. (1)-(8). For the case of labyrinth seal calculations, 
the flow information obtained often consists of a mass flux entering the 
labyrinth seai inlet chanvel, the entry temperature, Te, an entry 
pressure, Po, and the seal geometry. From the seal geometry an arbitrary 
reference length, 2, can be chosen (often che clearance gap height or 
the iale. channel height, 4.e., rotor to land distance is used). Given the 
inlet or free stream tempevature, Eq. (6) can be used to calculate the free 
stream viscosity, Wo. Since the mass flux, PwWeo, is “sown, Eq. (5) 
can be used to calcusate the Reynolds number per unit length, Re/ho. (In 
reality the mass flux will have to be multiplied by unity plus a blockage 
factor, b, to obtain the value of Poa. The dlockage factor is 
obtained from a knowledge of the inlet channel flow characteristics). 

Eqe (1) can then be usea to calculate the free stream density, pe, given 


the free stream pressure and temperature. Since the mass flux, 6,U,, 1s 


1+hb 
known the free stream velocity Us, can be calculated and nence the free 
stream Mach number, Mo, can be calculated. Finally at this point the 


flow conditions can be input via the wind tunnel option. 
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The MINT computer code has a very general coordinate capability, 


aso 


a 


i.e., calculations can be performed in general nonorthogonal coordinates. 
Both two-dimensional (ICOORD = 0) and axisymmetric (ICOGRD = 1) 
ncaorthogonal coordinate systems can be utilized. The former is primarily 


utilized to simulate laboratory test conditions while the latter is used to 


EL a Gs 


similate actual rotating hardware components. Geometric information can be 
calculated externally to the MINT computer code, in which case this data is 
read into the MINT code from unit NTAG = 21. In general, this capability 


is utilized for complex geometric configurations where it would be 


incoavenient to code such information into the HINT program. The MINT code 
has the capability of internally calculating the geometry for a wide 
variety of seal configurations. Basically, aay seal configuration 
consisting of knives and lands with either wall rotating that is 
constructed from a series of alternating nominally horizontal and vertical 
straight lines cac be -aualyzed. Thus, for example, rectangular, slanted 
and tapered kaives in both stepped seal and straight through seals can be 
analyzed (see Fig. 1). 

The proceas of constructing these ‘etick-like’ labyrinth seal 
configurations is as follews: (see Fig. 1). 

(1) Base lines (dotted lines of Fig. 1) are drawa on the bottom aad 
top and on the left-hand boundaries; 

(2) Grid points are next distributed on the bottom and top 


_ rae ES ae 


boundaries and on the left-hand boundary by use of the transformation 
technique of Oh (Ref, 2). Thies techuique allows selected grid points to be 


vs 
Pa Maa 
e 


conetrained to a specific physical location and thus the surfaces of tie 
rotor and knives can be specified; 


3, 
_ 


(3) Lines emanating from the grid points on the left-hand aide are 


cat TF. 
es Pil ee 3 


then drawn horizontally acrose to the right-hand side of the domain; aad 


Sone 
a. 


os 


(4) Grid pointe on th: cop and bottom boundaries are then connected 
with etraight lines thus completing the coordinate system. 

The method of distributing grid points on the boundaries fa described 
in detail in the next section. The basic procedure fa to concentrate a 
relatively large number of grid pointe in regions of the flow where large 
gtadients of the physical variablea are expected. For labyrinth seal 
configurations this usually seana in the vicinity cf boundary layere, ina 
the vicinity of the gap between the knives and the land snd in the 
vicinity of shock waves (if they exist). In general, {t is destrable to 
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start the computation domain far upstream of the first knife, i.e. on the 
order of four or five channel heights upstream of the knife. Since there 
will often exist a large recirculation zone downstream nf the last knife, 
the downstream region has to include at least the entire extent of this 
recirculation zone. The procedure normally employed is to use relatively 
large streamwise grid spacing at the upstream and downstream boundaries, 
and relatively fine streamwise resolution in the vicinity of the seal 
assembly. 

A wide variety of boundary conditions are available as input to the 
user of the MINT computer code. In addition, there are default values 
of boundary conditions which can be activated by identifying the boundaries 
with the FORTRAN variable IBOUND(KSURF,IDIR) where KSURF is the surface 
number with respect to the IDIRth direction (for discussion of convention 
used for KSURF and IDIR see section on Protocol for Grid Point and Boundary 
Identification). Values of IBOUND = 1,2,3 or 4 correspond to an inlet, a 
symmetry surface, an exit plane or a wall, respectively. Default values of 
the boundary conditions (FORTRAN variable IEQBU(KSURF,IDIR,IEQ)) associated 
with specific values of IBOUND are described in the NAMELIST input 
description section. Any default value of IEQBC can be overwritten by the 
inputing of a value of IEQBC in the NAMELIST input. 

One boundary condition that needs some explanation is the so-called 
‘two-layer model’ which is used as one of the default boundary conditions 
on an inlet. This boundary condition assumes that the flow in an inlet is 
divided into two regione or layers, a central inviscid core where the 
stagnation pressure is constanc., and an attached boundary layer(s6) where 
the preseure ia constant and the velocity profile has a prescribed foro. 
The prescribed form (but not the magnitude) of the boundary layer velocity 
profile must be supplied by the user. The MINT code presently has two 
methods of specifying the form of the boundary layer profile(s): 

(1) laminar profiles are generated in SUBROUTINE VKPOL where a 

von Karman~Polhausen (Ref. 3) profile {ts generated, and (2) turbulent 
profiles are generated in SUBROUTINE TURBP where a Maise-McDonald (Ref. 4) 
profile is generated. One feature of the two-layer model is that the mass 
flux will be determined from the solution of the governing equations, and 
not from an imposit’ n of the inlet profiles ot velocity and density. 

Since the steady state solution is obtained by time marching the 


governing equationa uni‘: the eolution ceases changing, an initial guesa 
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of all flow variables must be generated. Because of the wide variety of 
the geometric configurations considered, it would be impossible to develop 
a single technique for obtaining initial conditions. Rather, the approach 
taken in this study is to assume that initially the flow is at rest with 
all conditions equal to the stagnation values. Then the back pressure is 
lowered over a prescribed number of time steps until the desired downstream 
pressure is obtained. As the pressure pulse proceeds from the downstream 
boundary to the inlet, the flow is gradually accelerated or draw through 
the device until a steady state solution is obtained. The technique used 
is to gradually lower the back pressure, Pg, from the initial stagnation 


value, Po, to the desired back pressure, Pa, by use of the relationship 
Py = Pyli-f) + fF, (11) 


where 


- contr ma 0, nia(, st} (12) 


2 


and n is the present time step, nj is the time step when the back 
pressure, Pg, starts to change from ite initial value, P9, and nz, is 
the time step when Pg ie equal to the final pressure Pg. When 

n <a, f = o hence Pg = Po. When un > ng, f = 1 hence Pp & Pa. 

At intermediate values of n, f.e., ny <n ¢ 02, f varies as 


I+ cos (1) (13) 


2 


At values of n = ny and a = 02, se = 0. 

The choice of the numerical time step size fe critical in obtaining 
an economic converged solution to the governing equations. Since in this 
etudy only steady estate solutions were considered, there {6 no need to 
require that tranelent accuracy be maintained. One means that is used to 
rapidly obtain a steady state solution is to utilize a time step that 
varies with the location within the computational domain. The approach 
used here ie to use a sasll time atep in the region of the knives 


and larger time etepa in the regions both upstream and downstream 
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of the seal. This has the advantage of allowing the propagation of the 
pressure wave (caused by the lowering of the back pressure) to rapidly, 
i.e., in a few time steps, travel upstream. In the region of the knives, 
where large transverse and streamwise velocities and velocity gradients 
will occur, stability can be maintained with the use of the smaller time 
steps. In the region upstream of the knives, a larger time step is again 
employed as in this region the gradients can be expected to be small. 

The spatial time step variation, At, is controlled by use of the 


function ; 


At Aty | £(at-N+1+ 4 (at, -0 


(14) 
2 ~2, 
where 1 - cos4 wmax| 0, min(t, 7) 
f=t- 
eek 2 (15) 
z-2Z 
1 - cos 4 w max} O, min(t, 7) 
ine 4053 (16) 
20 2 


and where z,z},%2,23, and z4 are show in Fig. 2. At; is the ratio 

of the time step upstream of z) to the time step over the knife assembly, 
Sty, and Atg is the ratio of the time step downstream of Z, to the time 
step over the knife assembly. The variables f; and f2 serve a similar 
function to that of the previously described f in the back pressure 
formulation. 

The basic approach is to initially take aa large a time step as 
possible while still maintaining numerical stability. In this manner the 
basic features of the flow are rapidly developed and then the time step can 
be reduced to obtain the converged steady state solution. Typically the 
user must, for a given clase of geometric configurations, experiment with 
the initial time step. Usually the back pressure can be dropped over 
approximately 25 time steps, and with an additional 25 time steps the basic 
flow patterns will be developed. The important thing to observe during the 
initial time steps is the propagation of the pressure wave upstream. In 
addition, it fs desirable to watch the relative change (FORTRAN variable 
SSTST) of the dependent variables, and the values of the maximum residuals 
(FORTRAN variable RESMAX). Initially, both SSTST and RESMAX should be 
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relatively large. However, after the back pressure has been lowered to its 
final value both of these variables should start to decrease. By the 50th 
time step the streamwise velocity on the inlet plane should approach within 
15%-20% of its final value. The values of SSTST and RESMAX should be 
continuously monitored during a run as these are the parameters which are 
used to determine if the solution is converged. 

After the basic flow patterns have been developed, the large time 
step should be decreased to allow a converged solution to be obtained. 
As the MINT code is presently configured, the time step can only be changed 
during a restart. Hence, the overall strategy is to make an initial run at 
a large time step, to catalogue the restart files and then to wake 
subsequent runs at successively smaller time steps until a converged 
solution is obtained. The above guidelines have been successfully used in 
this study. However, it is important to stress that a certain amount of 
experience in time step selection is necessary, and several value 
judgements must be made in obtaining converged solutions. 

Solution of the governing. equations with a k-€ turbulence model 
requires the solution of two additional partial differential equations. 
The procedure used in this study is to first solve the flow equations of 
conservation of mass, momentum and energy with a mixing length turbulence 
model. The flow is then frozen in time and the k and € partial 
differential equations are solved to obtain a converged k and € field (for 
the frozen flow). Then the flow equations and the k and € equations are 
solved simultaneously until a converged solution is obtained. The above is 
done with the input FORTRAN variables IVISC and ISETKE. Initially IVISC is 
set equal to 3, the mixing length turbulence model option. A converged 
solution is obtained, and restart files catalogued. Next on a restart rua 
IVISC ig eet equal to -5, and ISETKE ie set equal to the last restart time 
step. An initial guess of k and € will then be obtained in SUBROUTINE 
SETDKE, and the k and € equations will be solved with a frozen flow field. 
Once a converged k-€ solution is obtained, a restart must be catalogued and 
anew run {initiated with IVISC set equal to 5 at which time the flow 
equations and the k and € equations will be solved simultaneously. 
It is posatble (but not recommended) that the intermediate step of 


solving the k-€ equations with frozen flow can be eliminated. 
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In this case the mixing length solution would be obtained and IVISC would 
be changed from 3 to 5 and the flow equations and the k and € equations 
would be solved simulataneously. The value of ISETKE would still have to 
be set to determine the time step when the transition from mixing length to 
k-€ turbulence model would occur. 

Solution of the governing equations with rotation requires that an 
addition swirl momentum equation be solved. The procedure used in this 
study was to first solve the governing equatiors with no rotetion and then 
gradually to increase the rotational speed of the appropriate boundaries 
from zero to the desired rotational speed. This is done by enforcing a 
boundary condition on the rotating surfaces. The rotational velocity on 


the boundaries is again obtained with the use of a cosine function of the 
form 


Vg = HONG (17) 


where r is the radius (at the desired point on the rotating surface), 
Veg is the final rotational speed and 


" { Sade n-n, 
cos ma ’ min{ Aa, |} (18) 
2 


where in this case n) fe the first time step where rotation starts and ng 


{cs 


is the time step when the final rotational velocity, Vege is achieved, 

It ts possible to save some computational time by not solving the 
rotational (swirl) momentum equation before time etep a, by setting 
ISWIRL = 0, and running until time step ny. At this point a restart must 


be catalogued and a new run initiated with ISWIRL = 1. This procedure is 
optional. 
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7.0 TRANSFORMATION FUNCTION 


The distribution of grid points on boundaries is accomplished by 
the use of a transformation technique developed by Oh and described in 
Ref. 2. If y and n designate the independent variables in physical and 
computational (grid point) space respectively, the transformation function 


f = on (19) 
can be integrated to yleld 
n 
y= f tndn + Yaug (20) 
Yenin 


A convenient transformation function is composed of a series of N 
complementary error functions of the fora 


t(y) = = By + th {erte[ 2 (»- % } 
. (21) 
ame (! + sign(a,)}} B, 


Tho qth couplementary error function is centered in computational space 
at location "py (which is referred to ae a pivot point), and a is 
the width in computational space ia which 90% of the grid size variation 
takes place. At computational space location “py. the values of £(n) 
wlll agsume a local waxiaus. Y is a convenient scaling conetant for ay 
{y = 1.163 x 2 = 2.326; erfc (1.163) = 0.10). In the limit as a, > 0, 
By is the yt etep height for pivot j, t.e., the difference in grid 
spacing on either side of Np 

Substitution of Bq. (21) into Eq. (20) and integration ylelde 


Y~ Yan = Agl(n-4,,,) + +i3{2 [ 4m) ~ 8heia)] 


Y 


~ [t+ atanta,l[a~ nau] A ba 
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where 
8, (yn) = = (9- Np) erfc [z (y- 1) 


_ b o-[r/a(n-ng, I (23) 
re e Pj 


The technique used in this study is to constrain the values of the 
physical coordinate, Your at specific values of the computational 
coordinate, Ne,> At interior points the Ncy,’s are referred to as 
interior cluster points. At the two end lines, there exist pairs of 
computational and physical points Nain, Yain» 29¢ "axe Ynax which 
are referred to as end cluster points. In this formulation (see Fig. 3), 
two pivot points are associated with each interior cluater point and one 
pivot point is asociated with each end cluster point (oain and 
\Daax’ Thus, if there are k interior cluster points, the total number 
of pivot points is N = 2k+2. To determine the functional relationship y(n) 
described by Eq. (22), requires (if the values of "py and ay are 
prescribed) that the values of the (N +1) B's in that equation be 
calculated. Both Eq. (21) and Eq. (22) are linear equations with respect 
to the 8's. By constraining the values of the physical coordinates at the 


interior and end cluster points, (k + 1) linear independent equations in 
_ the B'e are obtained. By further constraining the valuea of the slope, 
f(n), at the interfor and end cluster points, a further k + 2 Linear 
fadependent equations in the 6's are obtained and hence there are now 
L , 2k + 3 Linear tndependent equations for the (X + 1) =(2k + 3) B's. This 
system ig dolved for the 6's by standard Gauselan elimination techniques 


es ae Pe 9 


and hence all the constante of Eq. (22) are uniquely deteratned. 
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OPORTO 


There are several advantages to the use of the series of 


complementary error functions: 


(1) 


(2) 


(3) 


f(n) is positive, finite and non-zero, i.e., y will always 
increase with increasing n; 

f(n) is continuous and successively differentiable and 
integrable; and 

if the pivot points are spaced at a greater distance 

(in computational space) than a/2 from each other, the 
complementary error functions will have minimal interaction. 
Thus, the width and location of a complementary error function 
can be changed without affecting other complementary error 


functions. 


The strategy for using this technique is as follows: 


(1) 


(2) 


(3) 


(4) 


(5) 


Determine the location of the interior and end cluster points ia 
physical and computational space. The criteria used in 
determining these locations is usually based upon the physical 
processes that need to be resolved and the number of grid points 
needed to resolve those processes; 

The pivot points, Nps» and the corresponding bandwidth 
parameters, @ are then input. Usually the interior pivot points 
are located on either sides of the cluster points. If it is 
desired to have negligible interaction between the complementary 
error functions, the n*s and a's will be chosen such that no n 
will be located within a/2 grid pointe of each other; 

Slopea are then chosen at each intertor and end cluater point. 
Again these are usually determined from physical considerations; 
A test run is then made with a stand-alone version of SUBROUTINE 
OHGRID. Usually the initfal results will not be precisely to 
the liking of the user. In general this will be due to the 
choice of the slope constraints which can then be changed to 
obtain a better quality grid point distribution; and 

Once this interaction process has been coapleted with the slope 
constraints (because of item 3 in the advantages of this aethod 
that process ahould be rapid), slight modification of the values 
of the n‘s and a's should result in the desired grid point 
distribution. 
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8.0 PROTOCOL FOR GRID POINT AND BOUNDARY IDENTIFICATION 


Because of the complexity of the configurations considered in this 
study, a protocol had to be developed so that the user could easily set up 
the finite difference mesh required to solve the governing equations and 
apply the appropriate boundary conditions. The procedure utilizes the 
concept of a ‘grid array indicator' (FORTRAN variable IFSB (I, K), where I 
and K correspond to the x and z directions respectively) to type the grid 
points in the computational domain. Referring to Fig. 4 it can be seen 
that the idea is to identify the grid points according to the function they 
serve. The convention used to identify the grid point is as follows: 


IFSBC(I,K)=1 Indicates that a grid point is an interior fluid point. 


IFSB(I,K)=2 Indicates that a grid point is an interior solid point, and 
thus does not effect the calculation. 


IFSB(I,K)=3 Indicates that a grid point is a noncorner boundary point. 
IFSB(I,K)=4 Indicates that a grid point ia an inward corner point. 
IPSB(I,K)=5 Indicates that a grid point {6 a re~entrant corner point. 


Thus, for instance, grid points within the cross hatched areas of Fig. 4 
would all have values of I¥SB=2 while all the interior points would have 
values of IFSB@l. The boundary points (both on the cross hatched areas and 
on the outer boundartes) would have values of IFSB=3 except for the 


six fnward corner points (IFSB=#4), and the four re-entrant corner pointes 


(IPSB=5). The FORTRAN variable is dimensioned IFSB(71,101). Based upon 


the inforuation supplicd by the user, subroutine CORNER eete up the 
necessary patameters to control the ADL procedure. The grid array 
indicator is introduced through tIAMELIST S$READ6. 

In addition to the explanation of the grid array indicator, it {is 
convenient at Chis point to explain the surface number (KSURF) convention 
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and direction (IDIR) convention utilized throughout the MINT computer 
code, The correct data input for several variables, e.g., 
IBOUND(KSURF,IDIR), require a knowledge of this convention. The boundary 
surfaces are numbered according to a surface number relative to the 
computational direction. The convention is to allow values of IDIR=1 and 3 
to correspond to the x and z physical directions, respectively. Referring 
to Fig. 4, it can be seen that the surfaces (with respect to the 
x-direction) are numbered lx to 10x. This corresponds to values of KSURF 
ranging from 1 to 10 and a value of IDIR corresponding to the x-direction, 
i.e., IDIR=i. The analogous convention is also applied to the z-direction 
surfaces. Thus, the value of KSURF represents the surface number relative 
to the given direction, IDIR. 
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9.0 MINT INPUT 


Except for optional initial title cards, the entire MINT tnput is 
entered by means of the NAMELIST format. There are two primary advantages 
to the use of the NAMELIST format: (1) if the default values (defined by 
data statements) are acceptable, the user need not input that variable, and 
(2) the order (within a given NAMELIST) in which the input variables are 
entered is trrelevant. There are six NAMELIST input files in the MINT 
computer code, SREAD| through S$READ6. NAMELIST SREADI is used to enter 
testart information while NAMELISTs S$READ2 through SREAD6 respectively 
enter geometric data, flow conditions, time step and print control, 
boundary conditions and initial conditions. A detailed description of all 
the MINT computer code input will be given below. 
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9.1 


9.2 


Namelist or 
Variable Name 


SREAD1 


IREST 


IOTAPE 
INTAPE 


IOTAP1 


INTAPI 


$READ2 


NUMDX 


NUNDZ 


XGHIN( LDIR) 


XGHAX( IDIR) 


Description 


Restart Options 


Restart flag. 

IREST=0: New calculation 

IREST=1: Case is being run from a restart file 
Default value is 0. 


Output unit number for dependent variable array 
restart data 
Default value is 10. 


Input unit number for dependent variable array 
restart data 
Default value is 10. 


Output unit number for namelist restart data 
Default value is 20. 


Input unic number for namelist restart data 
Default value is 20. 


Geometric and Grid Options 


Number of interior grid pointe in the transverse 
direction (x direction). Total nuaber of 

pointe in this direction = NDUMKX + 2 

No default value. 


Number of interior grid points in the etreaawise 
direction (2 direction). Total nuaber of | - 
pointe tn this direction = NUMDZ + 2 

Ro default valwe, 


For IDIR = 1, disenstonless value of x-coordinate 
on bottos boundary. For IDIR © 3, dimensionless 
value of e~coordinate on left. boundary 

Default values are 3*C.0. 


For IDIR=1, dimensionless value of x-coordinte 

on top boundary. For IDIR=3, diaensionless value 
of e-coordinate on right boundary 

Default values are 3*1.0. 


IGEOM( IDIR) Geometry option 
IGEOM(IDIR)=1 - 


IGEOM(IDIR)=2 ~ 


IGEOM(IDIR)=3 —- 


IGEOM(IDIR)=4 - 


IGEOM( IDIR)=5 
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sentinel. 
Grid points equally distributed 
in direction IDIk. 


Grid points distributed in 
direction IDIR by method of Oh. 


Not used. 


Geometric data read from file 
NTAG=21. If this option is used, 
all values of IDIR must have 
IGEOM( IDIR)=4. 


Grid points distributed on lower 
and upper surfaces by method of Oh. 
IDIR=3 only. Coordinate system 
will consist of straight lines 
drawn from corresponding points 
on lower and upper surface. 
Remainder of coordinate system 
will consist of lines drawn 
parallel to lower (or upper 
surface) IGEOM(1) must then 
equal either 1 or 2. 


Default values are 3*l. 


LSHAPE Shape of computational domain 

L3HAPE*0 = Rectangular computational domain. 

LSHAPE=2~-4 =- Not used. 

LSHAPE=5 - Knives are input along lower 
surface. Location of the knife 
re~antrant corners are {nput 
via variables LCORNX and LCORNZ. 
Maximum number of corners ts 
NCRMAX=10, f.e., five knives. 

LSHAPE=6 - Computational domain is input 


through the use of the grid array 
indicator, IFSB (see NAMELIST 
SREAD6). 


Default value is 0. 
LCORNX( 1) Grid point number of x location of 1th cornet. 


No default values. 
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8 LCORNZ(T) Grid point number of z location of 1th corner. 


No default values. 


NCLUST( ISURF) Number of internal cluster points for SUBROUTINE 
OHGRID. When ISURF=1 points are distributed in 
x direction. Values of ISURF=2 and 3 
correspond to the distribution of grid points on 
the bottom and top surfaces respectively. Maximum 
number of internal cluster points equals 18 


a oo eee eee 


No default values. 


CLPX(I, ISURF) rth grid point number corresponding to the 
cluster points on surface ISURF. The first 
grid point on a surface must always be | while 
the last grid point on the transverse surface 
must be NDUMX + 2, and the last grid voint on 
the streamwise surface must be NDUMZ + 2. 

(1 =NCLUST(ISURF) + 2) < 20. 


No default values. 


CLPY(1, ISURF) The x or z coordinates (corresponding to ISURF#1 
and ISURF=2 and 3 respectively) which pair with 
the CLPX values. The first grid point on a surface 
must be the corresponding value of XGMIN(IDIR) and 
the last the corresponding XGMAX(IDIR) values. 


No default values. 


SLOPOH(I,1SURF) I¢h spacing constraint on surface ISURF. One 
value must be input for each CLPX(I, ISURF). 
However, if a negative value is input for the 
first and/or last values of CLPX(1,ISURF), 
no slope constraint will be enforced at these 
location(s). 


eee a oe 


No default values. 


Oe rs (Foe BF Fo 
e ” 


Ft. 


. a 


Bt PTAA IIMA 


33 


ee ee ee 


Se ee oe eit a Ca ghee Vtg Se TAA eA tate Sees 


Tons SA 4" Y a tite 


a rr 


ee) vate s oe ek a a a ie eee De Te ae Ce a ee Bee) a 
BO RI 1 0 OO INA Ba CN OO ES OT aT APOE ar 


ETAPOH(J , ISURF) 
; ALPHOH(J , ISURF) 
i 
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Jth jocation of a pivot point. At the first and 
last values of CLPX(I,ISURF), one ETAPOH(J,ISURF) 
must be input. Corresponding to each interim 
cluster point two ETAPOH's must be input. 

The values of ETAPOH are the grid point location 
of the centers of the series of complementary of 
error functions. (2*CLUSTC(CISURF) + 2) values of 
ETAPOH must be input for each value at ISURF. 
The rate of grid spacing will have the fastest 
variation at the pivot pointes. Pivot point(s) 
corresponding to the end points where values of 
SLOPOH are negative, are ignored, but must be 
input. 


J <¢ 38 
No default values. 


sth width apacing corresponding the location of a 
complementary error function at ETAPOH(J,ISURF). 
The width spacing is the number of grid points in 
which 902% of the grid spacing takes place around 
the pivot point. Again (2*CLUST(ISURF) + 2) values 
are required. Values of ALPHAH corresponding to 
the end points, where values of SLOPOH are 
negative, are ignored but must be input. 


J < 38 
No default values. 


Sentinel which determines if Cartesian or 
cylindrical polar formulation of governing 
system of partial differential equations are 
utilized. 


ICOORD=Q0 - Cartesian formulation. 


ICOORD=1 —- Cylindrical polar formulation. 


Default value is 0. 
In core ~ Out of core option. 
INOUT=0 - In core. 
INOUT=2 - Out of core option. 


Default value is 2. 


Flow Conditions 


System of unf{ts utilized. 
IUNITS=1 ~— MKS. 

IUNITS#2 - Engltah. 
Default value ie |. 
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CLENG Reference length (m or ft). 
No default. 

MINF Free stream Mach number. 
No default. 


TINF Free stream static temperature 
(°K or °R). Input required if REPL=0.0. 


No defauit. 
PINF Free gtream Stat{c pressure 
(nt/m” or lbf/ft~). Input required if 


REPL=0.0. If REPL#0.0, either PINF or 
PZERO must be input. 


No default. 
REPL Reynolds number per unit length (m7? or ft). . 
No default. 
PZERO Free stream stagnation pressure (Nt/m? or ine/ft*). 


No default. 


- 


IHSTAG Enthalpy option. 
IHSTAG*Q ~ Energy equation is formulated in terms 
of a static enthalpy. 


a 


= 
Fae a eS we 


~ 


IHSTAG=2 - Constant stagnation enthalpy is assumed. 
Default value is 0. 


ISWIRL Swirl equation option. 
ISWIRL=0 ~ No swirl equation is solved. 


| a tar tal l 


ged 2 


ISWIRL@1 - Swirl equation {8 solved - used for 


Re § cases with rotation. ICOORD must equal 1. 

. : Default value is 0. 

= ce! IVISC Viscosity model option. 

at IVISC#l - Conetant viecosity. 

N IVISC=2 - Laminar viscosity - Sutherland's law. 


IVISC#3 ~ Mixing length turbulent viscosity with 
Sutherland's laminar viscosity. 


IVISC=4 - k-fu turbulent viscosity with 
Sutherland's laminar viscosity. 
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9.4 


FK 


AVISC(IDIZ, IEQ) 


DIMIN 


DTMAX 


DT 


IDTADJ 


IVISC=5 - K-€ (Jones-Launder) turbulent viscosity 
with Sutherland's laminar viscosity 


IVISC=-5 - Solve K-€ equations only with frozen 
fluid dynamics. 


Default value is l. 

Wall roughness factor. 

Default value is 0.0. 

Artificial dissipation parameter I (See Eq. (1) of 
Appendix II) IEQ=1-5, 16 and 17 corresponding to 
the x,y,z momentum, continuity energy, k and ¢€ 


equations respectively 


Default values are all 0.5. 


Time Step and Print Control 

Number of time steps to be run. 

Default value is 0. 

Minimum nondimensional time step for this run. 

No default value. 

Maximum nondimensional time step for this run. 

No default value. 

Initial nondimensional time step. If DT is 
omitted on a restart, DT will be set to the value 
at termination of last run. 


No default value. 


Time step control. 
IDTADJ=0 - Constant DT is used for this run. 


IDTADJ«1 - Time step adjusted. If maximum relative 
change in any flow variable is less than 
0.04, DT is multiplied by 1.25 to a 
maximum of DIMAX. If maximum relative 
change in any flow variable is greater 
than 0.06, DT is divided by 1.25 to a 
minimum of DIMIN. 


IDTADJ@2 - Time step is cycled betwen DIMIN and 
OTMAX over NISTEP tiae steps. 


Default value is 0. 
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NISTEP 


NTREST 


IPRINT 


IVARPR(IV) 


Pekan * 
Snr ‘eee 8 
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Number of time steps used in cycling. 
Default value is 5. 

Restart is written every NIREST time steps. 
No default value. 


Complete flow field printout is provided every 
IPRINT time steps. 


No default values. 
Print control flag for variable IV. 
IVARPR(IV)=0 - Suppresses printout of variable. 


IVARPR(IV)#1 - Printout variable IV occurs. 


Iv IVARPR(IV) 
1 Transverse velocity - U 
2 Swirl velocity - V 
3 , Streamwise velocity - W 
4 Density ~ 9 
5 Enthalpy - h or hy 
16 Turbulence kinetic energy - k 
7 Dissipation of turbulence - € 
26 Static pressure ~ P 
27 Static temperature - T 
28 Effective viscoeity - vere 
29 Total temperature ~ Ty 
33 Mixing length - £, 
Kt) Dissipation - D:D 
35 Mach number - M 
36 Stagnation pressure - Py 


Default values are 1, 0, 1, 1280, 2%1, 8*0, 341, 
440, 1, 0, 21. 
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IGPRT(1V) Geometric variable printout option. 


eee Pe 


aoe 


IGPRT(IV)=0 - Suppresses printout of variable IV. 


IGPRT(IV)=1 - Printout of variable IV occurs 
on initial run. 


IGPRT(IV)=2 - Printout of variable IV occurs 
before first step of each run. 


Iv IGPRT(IV) 

ay) 

1 J = 

x} 

3 

oy” 

3 J 3x1 

ay? 

3 

9 1 ay” 

9x3 
13 r 
1? J 
21 x 
23 Z 


Default values are 1240, 341, 
O, 3*i, 0, 61. 


CSP 2 


IPLOT Plot option. 
IPLOT=0 <- No plot file written. 


eS a ty * 


IPLOT=999 = Plot file written on file TAPEL. 


C3 ele 


Default value is 0. 


3 

" ATIME1, ATIME2, ZDI, Input for spatial variation of time step - x 
Me 202, 203, ZD4 direction only. 

i See Eq. (14) 


Default value: are: 


ATIME1“0.0,  ATIME2=0.9 
201=10000.G ZD2=+10000.0. 
2D3=10000.G 2ZD4=-10000.0 


yan 


ISETKE Time etep when k-e equations sre first eolved 
No default value. 


ay Rrthy, AE OH 
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9.5 SREAD5 


IBOUND( ISURF , IDIR) 


TEQBC( ISURF, IDIR, IEQ) 
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Boundary Condition Input 


Type of boundary on surface number ISURF 
corresponding to direction IDIR. 


IBOUND( ISURF , IDIR)=1 


Inlet 
2 - Symmetry 
3 - Exit 
4 - Wall 


Maximum number of values in each IDIR 
direction is 50. 


No default values. 


Boundary condition on surface ISURF 
corresponding to direction IDIR associated 
with equation IEQ. IEQ has values of 
1,2,3,4,5,16 and 17 corresponding to x,y,z 
womentum equations, continuity, energy, k 
and € equations respectively. The 
defaulted values of IEQBC are determined by 
the values of IBOUND(ISURF,IDIR). 

For IBOUND (ISURF,IDIR)=1 the cross flow 
velocity componente are assumed to remain 
unchanged. The streanwise velocity 
component is determined from the two layer 
model. The pressure is extrapolated and 
the stagnation enthalpy remains constant. 
Values of k and € also remain unchanged. 
For IBOUND(ISURF,IDIR)=2 symmetry 
conditions are used, f.e., the first 
derivative of all variables except the 
normal velocity component are rero. The 
normal velocity component te zero. 

Por IBOUND(ISURF ,IDIR)=3 exit conditions 
are assumed. The second derivative of all 
velocity components, temperature, k and € 
are set to zero. The static presaure is 
preseribed. 

For IBOUND(ISURF,IDIR)=4 the no slip 
conditions are eet for the velocity 
components. The preesure condition is set 
by aesuaing thet the firet derivative of 
that variable is zero. The adfabatic 
condition fs assumed for the energy 
equation unless @ nonnegative value of 
TWALLCISURF,IDIR) is dnput, in which case 
the wall temperature {5 set to that value. 
Any default value IBQBC,(ISURP, LDIR, LEQ), 
can be overwritten by taputting a nonzero 
value of IEQBC. 
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TEQBC( ISURF , LDIR, LEQ) Boundary Conditions 


1 No change in function. 

2 First derivative equals 
Zere. 

3 Second derivative equals 
zero. 

4 Rough surface. 

19 Two-layer model momentum 
equation. 

50 Normal injection velocity 
specified. 

51 Normal masse flux specified. 

52 Rotational velocity 
specified. 

53 Rotational profile 
on inlet. 


z-direction grid point where injection begins. 
No default. 

z-direction grid point where injection ends. 
No default. 

Time etep when injection begins. 

No default. 

Time step when injection reaches final value. 
No default. 


Injection velocity (IEQBC=50) on injection mass 
flux (IEQBC#51). 


No default. 

Tie step when rotation begins. 

No default. 

Time step when rotation reaches final value. 


No default. 
40 
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i VROT Rotational speed (radians/sec). 

F No default. 

: NPRES1 Time step when back pressure begins to change. 

4 ‘ Default value is 0. 

e NPREP2 Time step when back pressure reaches final value. 
P 2 No default. 

; PAMB Ambient pressure. 

a No default. 

‘i TWALL(ISURF , IDIR) Dimensional wall temperature, °K or °R. 

Nae Default values are 150*0.0. 


gt 


CFI Skin friction coefficient on bottom surface of 
inlet. Needed for turbulent flow. 


No default. 


CF2 Skin friction coefficient on top surface of inlet. 
Needed for turbulent flow. 


Fi ETRE + fs, 


- 


No default. 


ow 
wr 


DELTAL Boundary layer thickness on bottom surface of 
inlet. 
No default. 

DELTA2 Boundary layer thickness on top surface of inlet. 


No default. 


TEMPUl Temperature of bottom surface, °K or °R. 
3 No default. 
TEMPW2 Temperature of top surface, °K or °R. 
No default. 
ORLTR Thicknese of rotational boundary layer on rotor. 


The rotational velocity (nendimensional) in the 
inlet will be a maximum of VROT*XGMIN(1) on the 
rotor and decrease by weans of a cosine function to 
zero at XGMIN(1) + DELTR. 


No default value. 
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Grid Indicator 


Indicator to type the grid points in the 
computational domain. I corresponds tg the x(y?) 
direction and K corresponds to the 2(y~) 
direction. Only used when LSHAPE=6. 


IFSB(I,K)=l - Interior fluids point. 


IFSB(I,K)=2 - Interior body point. 


IFSB(I,K)=3 - Boundary point. 
IFSB(I,K)=4 - Boundary point - re-entrant corner. 
IFSB(I,K)=5 - Boundary point - inward corner. 


Default values are 7171%*2. 
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10.0 SAMPLE INPUT 


Sample input for the MINT computer code is presented in Table 2. 
The case considered is that of a worn single knife seal (case 6a of 
Ref. 1). The first card of the input is an identifier for the input. The 
remainder of the input is entered through NAMELISTS $READ] through $READ6. 
This particular run is for a restart (IREST=1). However, except for the 
IREST#=1 card, the input for che initial run would be tdentical. NAMELIST 
SREAD2 indicates that the geometry utilized in this rua is obtained from 
file TAPE21 (IGEOM=3%4) and that a two-dimensional coordinate system is 
used (ICOOND=0). Values of NUMDX=29 and NUMDZ=99 indicate that 31 and 
101 grid points are used respectively in the x and z directions. 


SG ee Ee iis, tied is tee tat a. Rs 


The nondimensional dimensions of the physical domain are given by XGMIN and 
XGMAX as 12.0 and 300.0 in the x and 2 direction, respectively. <A value of 
LSHAPEeO indicates that no reentrant corners occur in the domain. NAMELIST 
S$READ3 gives the reference conditions through PORTRAN variables CLENG, 
MINF, TINF and REPL in English unite (IUNITS#1). A value of IHSTAG=0 
indicates that the energy equation te to be solved while a value of IVISC=5 
indicates that a k-€ turbulence model is used. NAMELIST SREAD4 indicstes 
that this case is to be run for 50 time steps (NT=50) from the last 
testart. This run is to be made at a constant time step of DT=0.05, 
DTMIN=0.05, DIMAX=0.05 and IDTADJ@O. Printout will occur each 25 tiae 
steps (IPRINT#=25), and variables u, w, k, €, p, T, Wage and M will be 
printed. ISETKE=650 indicates that this is the tiae step when the k-¢€ 
equations were first solved. NAMELIST $READS indicates that the upper and 
lower surfaces are balls and that the left hand side is an inlet and the 
right hand side is an exit (IBOUND=2%4, 480, 50*0, 1,3,48%0). The default 
boundary conditions are used. Values of PAMB=0.5, NPRESI=0 and NPRES2=50 
indicate that the back pressure was dropped from the total source pressure 


to @ value equal to 0.5 of the reference pressure over the first 50 steps 
of the run. CFl<CF2=5.0E-03, DELTAL=@DELTA2<=0.5 and TEMPWI=TEMPW2 are the 
paraeeters used for the boundary layer profile on the inlet surfaces 
(rotor and land, reepectively). NANELIST SREAD6 is empty as the LSHAPE=0 
ie used. However, the SREAD6 and SEND cards must still be input. 
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11.0 SAMPLE OUTPUT 


Sample output is presented in Table 3. The output consists of 
(1) the time step history of the run, and (2) the resulting output. 
Only sample portions of the output are presented as additional output would 
be repetitive and serve no purpose. 
The time step history output consists of pertinent output needed to 
monitor the run to determine the performance (convergence) of the run. It ‘ 
consists of the time step information (time step number, total CP time for 
this run, the physical elapsed time, the time step, the viscous stability 
criterion and the value of IDTADJ). The next two lines contain the maximum 
normalized change during a time step (SSTST) and its grid point location. 
The remainder is the maximum normalized change for all variables. The last 
portion of the time step history output contains the residual information. 
This residual information consists of the maximum average and peak 
residual for all equations and of the average and maximum residual, with 
its grid point location, for each equation. 
The second portion of the output consists of che various output as 
deterained by the FORTRAN input varible IVARPR. For economy only the 
Streanwise velocity, W-VEL, and the turbulence kinetic energy, TKE, are 
presented. 
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12.0 PLOT PROGRAM USER'S MANUAL 
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The MINT computer program does not have a plotting capability. 
Rather a plot file is written in SUBROUTINE WRPLOT and the user is free to 
interface with any desired plotting program. A brief description and 
user's manual for a plotting program, TDPLOTS is presenced below. Pregram 
TDPLOTS was specifically developed to interface with the MINT computer 
code, hence the plot file written in SUBROUTINE WRPLOT directly interfaces 
| with TDPLOTS. This code has the capability of producing plots of (1) the 

coordinate system, (2) vector plots of the velocity components, (3) contour 
" plots of up to 10 variables and (4) profile plots of up to 10 variables. 
TDPLOTS has a general three-dimensional capability, however in this. 
% description only the two dimensional option will be discussed. All 
plotting in program TDPLOTS is done in an x-y coordinate system where, 
by convention, x is the abcissa and y is the ordinate. In the discussion | 
of the input variables, that convention will be maintained. This progran 
has the capability of zooming in on any rectangular area of the physical 
domain. The input of FORTRAN variable control the zcom capability so that 
ic is not hardware liuited. The subroutines of prograa TDPLOTS are uritten 
so that they are not limited to a particular graphics library. This is 
¥ done by having all interface logic pase through one subroutine, PLTPKG, 
which aust be specialized for each lfbrary. At present four versions of 
a subroutine PLTPKG exist and can be used on the following graphics library 
4 Byetems: (1) CALCOMP, (2) IBM 3033 = NASA Lewis, (3) PLOTPSE and 
. (4) META. The toput for program TDPLOTS is input by asane of the NANELIST 
format. Three NAMELISTS SPLOT, SPLOT2 and SPRFPLT ate utilized. A 
deactipiton of the tnput associated with these NAMNELIST's will be given 
below. This will be followed by sample input. 
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i. 12.1 §PLOT Control Input 
4 RCASE Number of cages to te run. The number of cases is 


x defined aa the number of zooas required. In each 
‘ case fnforaatton must be tnput with NAMELIST $PLOT2 
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12.2 $PLOT2 


ICOORD 


IVECT 


VECM 


IKILLX(L) 


IKILLY(J) 


IPLOT(IV) 


Case Options 

Coordinate plot option. 

ICOORD=0 - No plot. 

ICOORD=1 - Plot. 

Nefault value is l. 

Velocity vector plot option. 

IVECT=0 — No plot. 

IVECT=1 - Plot. 

Default value is 0. 

Scaling for velocity vectors; length of tail of 
velocity vectors will be scaled by VECM. If 
VECM=0.0, scaling will be done on maximum value of 
magnitude of the velocity in the computational 
domain. 


Default.value is 0.0. 


Parameter to cull specific velocity vectors in the 
transverse direction at Ith grid point. 


IKILLX(I)=0 - Cull 

IKILLX(I)=1 - No cull. 

Default values 101*1. 

Parameter to cull specific velocity vectors in the 
normally streamwise direction at the Jth grid 
point. 

IKILLY(K)=0 - Cull. 

IKILLY(K)#1l - No cull. 

Default values are 101*1. 

Contour plot option 

IPLOT(IV)=0 - No plot. 


IPLOT(IV)=1l - Plot. 


Default values are 10*0. 
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THO(IV) 


HTS(I,1V) 


ISCL 
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Iv Contour Parameter 
1 U 
2 Vv 
3 W 
4 f°) 
5 P 
6 T 
7 M 
8 » (stream function) 
9 Cp 
10 Cp, 


Number of contour levels for variable IV. 
Maximum 10. 


Default values are 10*10. 


Sentinel to deteruine if contour levels are to be 
daput or calculated. 


THO(IV)=O - Contour levels calculated. Maximum and 
minimum levels of the contour variable are 
determined by scanning the computational domain. 
Plotted contour levels will then be equally spaced 
between the maximum and minimum at NHTS(IV) - 
levels, 


LHOCLV)=1 = Contour levels input from variable 
HTS(I,1V). 


Default values are 10*0, 

Contour levels; value of I goes from I-NHTS(IV). 
Default values are 100*0.0. 

Zoom parameter. 

ISCL=0 ~ No zoom. 

ISCL=1 - Zoom. 


Default value is 0. 
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XSCLI 


XSCL2 


YSCLI 


YSCL2_ 


TXB] 


IXB2 


IY) 


First x coordinate location for zoom (used only if 
ISCL=1). If XSCLI1 is positive, plot (coordinate, 
velocity or contour plot(s)) start at XMN+XSCLI on 
the x axis. Here XMN is the minimum value of x in 
the field. If XSCL1 is negative, plot will start 
at XMX-XSCLI where XMX is the maximum value of x in 
the field. 

Default value is ~1x105°, 

Last x coordinate location for zoom (used only if 
ISCL=1). If XSCL2 is positive, last x coordinate 
is XMN+XSCL2. If XSCL2 is negative, last x 
coordinate {8s XMX-XSCL2. 

Default value ts 1x10°°, 


First y coordinate location for zoom; analogous to 
XSCL1. 


Default value is -1x10°°, 


Last y coordinate location for zoom; analogous to 
XSCLI. 


Default value ts 1x105°, 

Sentinel to draw lower boundary of plot domain. 
IXB1=0 - No draw. 

IxBle1 - Draw. 

Default value {fe i. 


Sentinel to draw top boundary of olot 
domain. 


IXB2=0 - No draw. 

IXB2=]1 = Draw. 

Default value is 1. 

Sentinel to draw left hand boundary of plot domain. 
IYB1@0 - No draw. 

IYBi=1 - Draw. 


Default value is 1. 
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TYB2 Sentinel to draw right hand boundary of plot 
domain. 


IYB2=0 - No draw. 

IYB2=1 ~ Draw. 

Default value is l. 
NPRFIG Number of profile plots. 


If NPRFIG=0, NAMELIST PRFPLT input is omitted. 
Maximum number of figures is 10. 


Default value is 0. 


12.3 $PRFPLT Profile Plot Input 
NPRF( IP) Number of profile plots on IPth figure. 


Default values are !0*1. 
NVRPL(IV,IP) Variables to be plotted on IP‘ profile plot. 


Default values are 100*0. 


Variable NVRPL(IV, IP) 


U 1 
v 2 
W | 3 
p 4 
P > 
: 6 
M 7 
v 8 
Cp 9 
py 10 


CRDIC(I, IP) CRDI(L,IP) ~- x coordinate of the first profile 
point of the rPth plot. 


CRDI(2,I1P) ~ x coordinate of the last profile point 
of the IPth plot. 


No default values. 
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CRD2(1L,IP) CRD2(1,IP) - y coordinate of the first profile 
point of the IPth plot. 


CRD2(2,IP) - y coordinate of the last profile point 
of the rPth plot. 


No default values. 

NPTS Number of equally spaced points used on straight 
line joining first and last profile points to 
obtain the interpolated values for the profile 
plot(s). Maximum value of 101. 

Default value is 50. 

XOFSET(IP) Constant to be added to profile plot independent 
variable. If XOFSET™0 independent variable will 
vary from 0 to the length of the line joining the 
first pofile point (CRDI(1,IP), CRD2(1,IP)) to the 
last profile point (CRD(2,IP), CRD2(2,IP)). 
Default values are 10*0.0. 

DSCALE(IP) Scale factor used to multiply dependent variables 
for IPth plot. Often used to plot dimensional 
values. 

Gefault values are 10*1,0. 

XLEN(IP) Length (in) of x axis. 
Default value is 10*6.0. 

YLEN(IP) Length (in) of y axis. 
Default value is 10%6.0. 

ITURN(LP) Flag used to rotate 1Pth profile plot. 
ITURN(IP)=O = Dependent variable along x-axis. 
ITURN(IP)=1 ~ Depandent variable along y-axis. 


Default values are 10*0. 
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PIMIN(IP), PIMAX(IP) 


PDMIN(IP), PDMAX(IP) 


Extents of independent variable axis for rpth 
plot. Note that the inch increment on the x-axis 
will be equal to (PIMAX(IP) - PIMIN(IP))/XLEN(IP) 
of independent variable. 


Default is automatically calculated using the 
independent variable values to be plotted. 


Extents of dependent variable axis for IPpt® plot. 
Note that the inch increment on the y-axis will be 
equal to (PDMAX(IP)- PDMIN(CIP))/YLENCIP) of 
independent variable. 


Default is automatically calculated using the 
dependent variable values to be plotted. 
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12.4 Sample Input 


Sample TDPLOTS input is presented in Table 4. The NAMELIST $PLOT 
input NCASE=2 denotes that two cases are to be run thus requiring that 
NAMELIST $PLOT2 must be input twice. The first NAMELIST $PLOT2 uses the 
zoom option (ISCL#=1), and will make plots between x-values of 12.0 to 18.0 
and y-values of 0.0 to 10.0. Since the default value of ICOORD is 1, a 
coordinate plot will be generated. The values of IPLOT will cause stendard 
contours (10 contour levels between the minimum and maximum field values) 
to be plotted for Mach number and stream function. No vector or profile 
plots are associated with the first SPLOT input, i.e., default values of 
IVECT and NPRFIG are both zero. The second NAMELIST $IPLOT2 input will 
retain the previously defined ISCL=1] value, and will plot contours of Mach 
number and stream function between x-values of 18.0 to 24.0 and y-values of 
0.0 to 10.0 (the values from the first S$PLOT2 input are retained). In 
addition, vector plots will be made cver the same domain (IVECT#=1). 
However, the coordinate plot will not be made (ICOURD=0). Two profile 
plots will be made (NPRFIG=2) necessitating S$PRFPLT input. Both plots will 
contain two profiles (NPRF(1)=2*2) of u and w velocity components 
(NVRPL(1,1)#1,3,8%0,1,3,8%0) at a constant x station of 20.0 with y-values 
varying between 0.0 and 10.0. 
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LIST OF SYMBOLS 


Symbols 

B Constant from quadratic equation (Eq. 9) 
; b Blockage factor 

Cp Specific heat at constant pressure 

cy Constant in Sutherland's viscosity law 

c2 Constant in Sutherland's viscosity law 

f Cosine factor or transformation function 

h . Enthal py 

£ Length 

M Mach number 

n Time step number 

P Pressure : 

R Universal gas constant 

Re Reynolds number 

r Radius 

tT Temperature 

t Time 

U Velocity 

Ve Rotational velocity 

y Physical distance 

z£ Streamwise distance 


Greek Syabola 


@ Computational width factor 

] Step height for pivot 

Y Ratio of apecific heats or sealing constant 
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Subscripts 


“ac, 


a 
B 
f 
j 
y o 
. P 
*% 
1 
x 1 
.. 2 
, : 
x 4 
Mt] 
i max 
% min 
: - 


SBT ISL IS 5 '4 ca 


Over 


Change 


Computational distance 


Grouping defined Eq. (23) 


Viscosity 


Density 


Ambient 

Back or knife assembly 
Final 

th 

Stagnation; first 
Pivot 

First 

Second 

Third 

Fourth 

Maxinun 


Miniaun 


Free stream or reference condition 
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Pigure 1 - Typical Coordinate System Construction. 
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Figure 2 - Time Step Selection. 
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<a 
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. *kIRKWORN SINGLE KNIFE SEAL*¥*** 
mt SREADI 
in IREST = 1, 
ott SEND 
“ $READ2 


NUMDX = 29, NUMDZ = 99, 


IGEOM =» 3%*4, ICOORD = 0, 
XGMIN = 0.0, 0.0, 0.0, 
XGMAX = 12.0, J.0, 300.0, 
LSHAPE = 0, 
. SEND 
$READ3 
IUNITS = 1, 
CLENG = 2.54E-03, 
REPL = 2.4309269E+06, 
MINF = 5,4357E-02, 
PZERO = 2.004087E+05, 
' IRSTAG = 0, 
IVIsc = 5, 
SEND 
SREAD4 
NT = 50, 
IPRINT = 25, 


PAMB = 0.5, NPRES! = 0, NPRES2 = 50, 
ISETKE #650, 

SEND 

SREADS 

IBOUND = 2*4,48%0, S0*0, 1,3,48%0, 


TEMPWL = 294.44, TEMPW2 = 294.44, 
SEND 

SREAD6 

SEND 


« 


RD cee 


, ieee 
*e "e 


Pie 
a, Eo 


Table 2 - Sample MINT Input 


IDTADJ = 1, DT = 0.05, DIMIN = 0.05, DIMAX =0.05, 
IVARPR = 1,0, 1, 12*0,2*1,8%0, 31,640, 21, 


DELTAL = 0.5, DELTA2 = 0.5, CF] © 5.0E-03, CF2 « 5.0E-03, 


- Sample MINT Output 
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